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Abstract 

We study bottom-quark fragmentation in e~^e~ annihilation, top and Higgs 
decay H — > bb, using Monte Carlo event generators, as well as calculations, 
based on the formalism of perturbative fragmentation functions, which 
resum soft- and collinear-radiation effects in the next-to-leading logarithmic 
approximation. We consider the PYTHIA and HERWIG generators, and 
implement matrix-element corrections to the parton shower simulation of 
the H bb process in HERWIG. We tune the Kartvelishvili, string and 
cluster models to 5-hadron data from LEP and SLD, and present results in 
both xb and moment spaces. The S-hadron spectra yielded by HERWIG, 
PYTHIA and resummed calculations show small discrepancies, which are 
due to the different approaches and models employed and to the quality of 
the fits to the e~^e~ data. 



1 Introduction 



Heavy-flavour phenomenology is currently one of the most lively fields of investigation, 
in both experimental and theoretical particle physics. In this paper, we consider B- 
hadron production in e+e~ annihilation {e'^e~ — > bb), top decay {t — > bW) and the 
decay of the Standard Model Higgs boson H — > bb. In fact, data on 6-flavoured hadron 
production from SLD [1] and LEP [2,3] experiments have been available for a few years, 
and they can be used to predict 6-quark fragmentation in other processes. Bottom 
fragmentation in top decay is indeed one of the main sources of uncertainties on the 
top-mass reconstruction at the Tevatron [4] and the LHC. In particular, the analysis in 
Ref . [5] , where the top-quark mass is determined from the invariant mass of three leptons 
(one lepton from the W decay and two leptons from a J/ip decay), depends on a 

proper description of the fragmentation of h quarks into B hadrons which decay further 
to J/'^'s. As for H — > bb, the detection of this process is impossible in the channel 
gg{QQ) — > — > at hadron coUiders, because of the large irreducible bb background. 
The solution to this problem is to look for H bb in associated production channels. 
The most relevant search channels are Higgs production in association with a vector 
boson at the Tevatron [6], and in association with a W boson [7] or with ti pairs [8] at 
the LHC. In such channels, the H ^ bb decay mode can be probed for a Higgs mass 
up to niH < 135 GeV. For heavier Higgs bosons, the decay to bb might be accessible in 
Two- Higgs Doublet Models, with enhanced Higgs couplings to bottom quarks [9]. In all 
cases, a good understanding of bottom fragmentation in Higgs decays is important for 
the Higgs mass reconstruction and the identification of 6-jets, which both depend on the 
properties of the Higgs decay products. 

In order to perform accurate searches and measurements, the use of precise QCD 
calculations will be fundamental. In particular, fixed-order calculations are reliable 
enough to predict total cross sections or widths, but differential distributions exhibit 
terms, corresponding to coUinear or soft parton radiation, that need to be summed to 
all orders to obtain a meaningful result. An example is the large term ~ ln((5^/m^), 
where Q is the hard scale of the process and m the heavy-quark mass, which appears in 
the energy distribution of a heavy quark at next-to-leading order (NLO). 

A possible approach to study 6-quark production in the considered processes is based 
on the perturbative fragmentation formalism [10]. The NLO heavy-quark energy spec- 
trum is expressed as the convolution of a coefficient function, describing the emission of a 
massless parton, and a perturbative fragmentation function D{mh,iiF), associated with 
the transition of a massless parton into a massive quark. The dependence of D{mb, /ip) 
on the factorization scale hf is determined by using the Dokshitzer-Gribov-Lipatov- 
Altarelli-Parisi (DGLAP) evolution equations [11, 12], once an initial condition at a 
scale yUoF is given. The initial condition of the perturbative fragmentation function, 
first computed in [10], has been proved to be process- independent in [13]. Solving the 
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DGLAP evolution equations we can resum the large ln{Q'^/m'^) (coUinear resummation) . 

Moreover, heavy-quark energy distributions present terms that are enhanced when 
the energy fraction x approaches 1. Such contributions can be resummed following the 
general method of [14, 15] (soft or threshold resummation). Soft and coUinear resumma- 
tions in the next-to-leading logarithmic approximation (NLL) have been implemented 
for b production in e~^e~ annihilation [10,13], top-quark [16,17] and Higgs [18] decays, 
in the framework of perturbative fragmentation functions. 

An alternative approach to address heavy-quark production and resum terms corre- 
sponding to soft and coUiner radiation consists in using Monte Carlo event generators. 
Standard Monte Carlo programs, such as HERWIG [19] or PYTHIA [20,21], simulate 
multiple emissions in the soft or collincar approximation, and arc provided with matrix- 
element corrections [22, 23] to describe hard or large-angle radiation. 

In order to describe hadronic spectra, perturbative calculations and Monte Carlo 
parton showers need to be supplemented by models describing the hadronization. Cal- 
culations based on the fragmentation-function approach typically predict hadron-level 
energy distributions convoluting the partonic spectra with a non-perturbative fragmen- 
tation function, such as the Kartvclishvili [24] or the Peterson [25] models. These models 
contain parameters that need to be fitted to experimental data. More recently [26], it 
was suggested that one can fit directly the A^-moments of heavy-hadron cross sections 
and extract the moments of the non-perturbative fragmentation function, with no need 
to assume any functional form in x-space. Likewise, Monte Carlo event generators sim- 
ulate the hadronization using appropriate non-perturbative models, such as the cluster 
model [27] and the string model [28], for HERWIG and PYTHIA, respectively. 

The outline of this paper is the following. In section 2 we review coUinear and soft 
resummations in the framework of perturbative fragmentation functions. In section 3 
we discuss Monte Carlo parton shower algorithms. Section 4 describes matrix-element 
corrections to parton-shower simulation and, in particular, the implementation of such 
corrections to if —> 66 processes in HERWIG. In sections 5 and 6 we shall present 
our results on S-hadron production in e'^e~ annihilation, top and Higgs decays, in x 
and moment space respectively, after tuning the hadronization models to e'^e~ data 
from LEP and SLD. Section 7 summarizes the main results and gives some concluding 
remarks. 



2 Perturbative fragmentation functions 

In this section we shall review the main points of the perturbative fragmentation ap- 
proach and of the implementation of coUinear and soft resummations. We shall consider 
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b production at next-to-leading order in Z, top or Higgs decays, i.e. 

PiQ) ^ bip,)bipi) {g{p,)) , (1) 

with P = Z or H, and 

m^b{p,)W{pw){9{Pg)). (2) 
The energy spectrum of a massive b quark is given by the following general result: 



1 dr 



5(1 -Xb) + 



Q2 

Pqq{Xb)\n — + A{Xb) 



Tq dxb ^ ' 2% 

In Eq. ©, 3^6 is the normalized 6-quark energy fraction: 

1 2pb-Q 



Xb 



w 



(3) 



(4) 



where w = in Z and H decay, and w = m^/ml in top decay, since the W mass is 
not negligible with respect to the top mass. Furthermore, in Eq. (jS)) Fq is the width of 
the Born process, /i is the renormalization scale, p > 1, Pqq{xb) is the Altarelli-Parisi 
splitting function: 



Pqq{Xb) = Cf 



l-Xb, 



(5) 



Function A{xb) is process-dependent and does not contain the b mass. 



The large logarithm ln(Q^/m^), which appears in Eq. (jS)), can be resummed following 
the approach of perturbative fragmentation functions [10]. The b spectrum is expressed 
as the convolution of a massless coefficient function and a perturbative fragmentation 
function D{mb, fip), associated with the transition of a massless parton into a heavy 
quark: 



Fo dxb 
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(6) 



In Eq. dUj), dVi/dz is the differential width for the production of a massless parton i in 
any of the considered processes, after subtraction of the coUinear singularity in the MS 
factorization scheme. Throughout this paper, and for all the considered processes, we 
shall neglect secondary 6-quark production via g ^ bb splitting, which means that, in 
Eq. i = b and expresses the fragmentation of a massless b into a massive b. 
The NLO MS coefficient functions were calculated in [10,16,18] for Z, t and H decays 
respectively. 



The perturbative fragmentation function follows the DGLAP evolution equations 
[11, 12]. Its value at a given scale hf can be obtained once an initial condition is given. 
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In [10] the initial condition D™\xb, fioF, was calculated and its process-independence 
was established on more general grounds in [13]. It is given at NLO by: 



asii^DCp \ l +xl 
2% I — Xb 




21n(l-Xb) - 1 



J + 



(7) 



As discussed in [10], solving the DGLAP equations for an evolution from fiQp to 
fiF, with a NLO kernel, allows one to resum leading (LL) ln"(/i|.//iQ^) and next- 
to-leading (NLL) ln"~^(/i|'//^Qp) logarithms (collinear resummation) . The explicit 
expression for the solution of the DGLAP equations can be found, for instance, in 
Ref. [10]. Setting //qf — f^b and /i^ ~ Q, one resums the large \n{Q'^/ml) appearing in 
the massive spectrum Q. 

The NNLO initial condition of the perturbative fragmentation function was calcu- 
lated in Refs. [29,30]; however, for the NNLO result to be applicable, one would also 
need three-loop time-like splitting functions, which are currently unknown. 

Furthermore, both the initial condition (|7j) and the coefficient functions in [10,16,18] 
present terms, ~ 1/(1 — 0;;,)+ and ~ [ln(l— a;b)/(l — 0;^)] + , which become large for Xb — > 1. 
The large-Xfe limit corresponds to soft-gluon radiation; in Mellin moment space, such 
terms correspond, at 0{as), to single (~ as\nN) and double (~ a^ln^A^) logarithms 
of the Mellin variable N. Soft resummation in the initial condition of the perturbative 
fragmentation function is process-independent and was performed in [13] in the NLL 
approximation. Large-x^ resummation in the coefficient functions was implemented 
in [13] for e^e~ annihilation, in [17] for top decay, in [18] for H ^ bb processes. In A^- 
space and in the NLL approximation, terms ~ ln"+^ AT (LL) and ~ agin" AT (NLL) 
are kept in the Sudakov exponent. 

The soft and collinear resummations in [13, 17, 18] are matched to the exact NLO 
results, in order to give a reliable prediction over the full Xb (N) range. This implies 
that the total widths are NLO as well. 

It was shown in Refs. [13, 17, 18] that, after soft and collinear logarithms are re- 
summed, the 6-quark spectra present very little dependence on factorization and renor- 
malization scales, which is equivalent to saying that the theoretical uncertainty is re- 
duced. Hereafter, we shall set the scales appearing in Eqs. (0) and ((Tj) to: fi = fip = Q, 
yUo = /^OF = mb. 
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3 Parton shower algorithms 



In this section we discuss parton shower algorithms, as implemented in Monte Carlo 
event generators, such as HERWIG and PYTHIA. These algorithms rely on the uni- 
versality of the elementary branching probability in the soft or coUinear approximation. 
Referring, for simplicity, to final-state radiation, the probability of emission of a soft or 
collinear parton reads: 



2t, ' ' AsiQ^Ql) 

In dHJ, P{z) is the tree-level Altarelli-Parisi splitting function ^, z is the energy fraction 
of the emitted parton with respect to the emitter, is the ordering variable of the 
shower. In HERWIG [19], is an energy- weighted angle ^, which corresponds to 
angular ordering in the soft limit [31]. In PYTHIA [20], is the momentum squared 
of the radiating parton, with an option to veto branchings that do not fulfil the angular 
ordering prescription. Moreover, the latest version of it, PYTHIA 6.3 [21], offers as 
an alternative the possibility to order final-state showers according to the transverse 
momentum of the emitted parton with respect to the emitter's direction, along the lines 
of [32]. For most of the results that we shall present, we shall use PYTHIA 6.220, with 
the option to reject non-angular-ordered showers turned on, and HERWIG 6.506 [33]. 

In (jS)) As{Ql, Ql) is the Sudakov form factor, expressing the probability of evolution 
from Qf to Q2 with no resolvable emission. In particular, the ratio of form factors in 
(jHl) represents the probability that the considered emission is the first, i.e. there is no 
emission between and Q^ax) where Q^ax is set by the hard-scattering process. Ql is 
instead the value of at which the shower evolution is terminated. In diagrammatic 
terms, the Sudakov form factor sums up all virtual and unresolved real emissions to all 
orders. 

For multiparton radiation, iterating the branching probability (jH)) is equivalent to 
performing the resummation of soft- and collinear-enhanced radiation. As discussed, 
for example, in [34] in the framework of the HERWIG event generator, parton shower 
algorithms resum leading logarithms in the Sudakov exponent, and include a class of 
subleading NLLs as well. 

Standard Monte Carlo event generators yield the leading-order total cross section 
for all implemented processes. The more recent 'Monte Carlo at next-to-leading order' 



-'^ Unlike Eq. which includes virtual corrections as well, in Eq. ((HJ P{z) only accounts for real, 
small-angle parton radiation. We have, for example, Pqq{z) — Cf(1 + 2:^)/(l — z), without the plus 
prescription appearing in Eq. 0. 

^In the HERWIG showering frame, ~ E'^{1~ cos 9), where E is the energy of the splitting parton 
and 9 is the emission angle [31]. 
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(MC@NLO) [35] program implements instead both real and virtual corrections to the 
hard-scattering process, in such a way that predicted observables, including total cross 
sections, are correct to NLO accuracy. 

While resummed calculations typically do not implement the widths of the decaying 
particles, whose masses are fixed quantities in the calculations, Monte Carlo generators 
include the finite widths of Z, W, Higgs bosons and top quarks. Some discussion on 
width effects will be presented in section 6, when we shall investigate the decay of a 
Higgs boson with a mass of 500 GeV. 

Moreover, for Higgs or top quark events, both HERWIG and PYTHIA neglect in- 
terference effects between production and decay phases. In this approximation, the 
differential decay width, normalized to the total width, is equal to the normalized dif- 
ferential cross section, independently of the production process: 

1 dT 1 da 
r dxb (J dxb 

Neglecting interference is a reliable approximation, as long as the energies of the radiated 
gluons are much larger than the top or Higgs widths [36]. This is indeed the case 
whenever the experimental analyses set cuts, e.g. on the transverse momenta of the 
reconstructed jets in top-quark events, of the order of 10 GeV or larger. Although not 
really essential in view of Eq. (j^I), in the following we shall nonetheless run HERWIG 
and PYTHIA for top and Higgs production at the LHC, i.e. for pp collisions at a 
centre-of-mass energy ^/s = 14 TeV. 



4 Matrix-element corrections to parton showers 

The algorithm briefly discussed in the previous section is able to simulate soft or collinear 
radiation; for hard and large-angle emission, the exact matrix element must be used. In 
fact, Monte Carlo parton showers are supplemented by matrix-element corrections. 

PYTHIA uses the soft / collinear approximation in all the physical phase space and the 
exact tree-level matrix element corrects the flrst emission [23,37]. In particular, PYTHIA 
6.220 contains matrix-element corrections to all processes that we shall investigate. 

The standard HERWIG algorithm entirely suppresses radiation in the so-called 'dead 
zone' of the phase space, corresponding to hard and large-angle emission. The exact 
matrix element is used to populate the dead zone (hard correction) and to correct the 
shower every time an emission is capable of being the 'hardest so far' (soft correction), i.e. 
the emitted parton has the largest transverse momentum with respect to the emitting 
one [22]. HERWIG 6.506 [33] contains matrix-element corrections to e+e~ annihilation 
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[38], Deep Inelastic Scattering [39], top decay [40] and vector-boson production [41]. 
The corrections to Higgs hadropro duct ion {gg/qq ~^ H), whose inclusion is in progress, 
are discussed in [42]. In this section, we shall briefly discuss the implementation of 
matrix-element corrections to parton showers in if —> 66 decays, not yet present in 
HERWIG. 

As the Higgs is a colourless particle, the corrections to H ^ bb are a straightforward 
extension of the ones to Z —>■ bb processes [38]. The limits of the phase-space region 
populated by HERWIG will be exactly the ones computed in [38], with ttih replacing 
rriz- The total and HERWIG's phase spaces are plotted in Fig. HI in terms of Xb and x^. 
Collinear, corresponding to either Xb = 1 or xi = 1, and soft emissions, i.e. the point 
Xb Xb 1 , are within the HERWIG phase space; no radiation is allowed in the dead 
zone. We also note an overlapping region, where radiation may come from either b or b. 



1^ 
X 




0.2 0.4 0.6 0.8 1 



Figure 1: Total (solid) and HERWIG (dashes) phase spaces for gluon radiation in H 
bb processes. The dead zone is empty in the standard algorithm. 



To implement matrix-element corrections, we need to compute the double-differential 
width of H ^ bbg processes, corresponding to the Feynman diagrams in Fig. in terms 
of the b and b energy fractions. It reads: 

1 d'^T _ asCp xl + xl + 2xbX-b - 2xb -2x^ + 2 
Fq dxb dxb 2n (1 — Xb){l — xj) 

In Eq. (|Tn|) the scale of as is set to the gluon transverse momentum with respect to 
the direction of the emitting quark: this choice is suitable to sum up a class of next- 
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to- leading logarithmic corrections [15]. When applying matrix-element corrections, we 
generate events in the dead zone according to Eq. (fTIH) . and use the exact result to 
correct the cascade in the already-populated region, any time an emission is the hardest 
so far. We have estimated that, after matrix-element corrections are applied, for a Higgs 
mass rriH = 120 GeV, about 4% of events have a gluon emission in the dead zone. For 
niH = 500 GeV, this fraction becomes about 3%, since, for a larger mass value, even 
the gluon transverse momentum is on average larger and therefore cesiq^) smaller. 




Figure 2: Feynman diagrams for H bbg. 




Figure 3: 6-quark energy spectrum in Higgs decay, according to HERWIG with (solid) 
and without (dotted) matrix-element corrections, on linear (a) and logarithmic (b) 
scales, for tjih = 120 GeV. 

Figures OH^ show the effect of matrix-element corrections on the 6-quark energy 
fraction, for mn = 120 and 500 GeV, the values that we will consider in the following. 
As expected, more events are generated at small Xb (corresponding to the dead zone of 
the standard algorithm) through the exact amplitude; this enhancement is compensated 
by having less events at large Xb- In fact, even after matrix-element corrections are 
applied, HERWIG, as well as PYTHIA, gives by default the total leading-order cross 
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Figure 4: As in Fig. but for a Higgs mass nin = 500 GeV. 

sections (widths). We also observe that the Xh spectra present a sharp peak for Xh ^ I, 
which corresponds to events where no gluon is emitted by the bb pair {xb = xi = 1). The 
fraction of events with Xf, = 1 depends on the shower cutoff Qo in Eq. (jH)), a user-defined 
parameter of HERWIG [19]. 



5 5-hadron spectrum in x^-space 

In this section we wish to present results on the 5-hadron spectrum in e~^e~ annihila- 
tion, Higgs and top decay. Our tools to describe 6-quark production will be, as discussed 
above, NLL-resummed calculations, HERWIG and PYTHIA. Given the different per- 
turbative accuracies, treatment of width effects and hadronization models which are 
implemented, we do not expect that NLL computations and Monte Carlo event gen- 
erators should necessarily agree. However, we find it very interesting to perform this 
comparison, since they are tools that are used for the analyses of 6-quark fragmentation, 
and it is useful to estimate the quality of the agreement between the predictions they 
make. Since the description of perturbative b production is different, in order for our 
comparison to be consistent, we shall have to tune independently HERWIG, PYTHIA 
and the hadronization model used in the framework of resummed calculations to the 
same data set. Such a consistency is crucial if we wish to extract non-perturbative in- 
formation from e^e~ annihilation and use it to predict 5-hadron production in other 
processes [26,50]. 

In order to describe _B-hadron spectra in the considered processes, we use experi- 
mental data on B production in e~^e~ annihilation, and rely on the universality of the 
hadronization mechanism. We shall account for data from SLD [1] and from the LEP 
experiments ALEPH [2] and OPAL [3] ^. Our results will be expressed in terms of the 

•^Contrary to Refs. [16-18], here we consider the OPAL data too. 
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B normalized energy fraction xb, which is the hadron- level counterpart of Eq. 



1 2pB-Q 

Xb = z -7^, 11 

where is the B four-momentum. As pointed out in section 2, when we use a resummed 
calculation that neglects powers of mf/Q^, w = m^r/rnl for top decay, while if = for 
the other processes. Since Monte Carlo event generators include some vn^jQ^ effects, 
w = m^/m^ —ml/m^ for top decays in HERWIG and PYTHIA. 

ALEPH and OPAL reconstructed only weakly-decaying B mesons, while the SLD 
sample contains some B baryons too; however, since the fraction of baryons is pretty 
small, we can safely fit all data together and, when running HERWIG or PYTHIA, we 
shall account for both mesons and baryons. As the resummed calculation yields a NLO 
total cross section (width), and the considered Monte Carlo programs only the LO one, 
we shall study normalized distributions, such as l/F dT/dxB, everywhere. 

As far as the NLL resummed calculation is concerned, up to power corrections, one 
writes the hadron-level spectrum as the convolution of the parton-level one, determined 
as discussed in Section 2, with a non-perturbative fragmentation function: 

i^(x.,Q,m.) = i r ^f^(.,g,m,)Z^- (^). (12) 
L axB 1 z az \ z J 

We describe the hadronization using the Kartvelishvili non-perturbative fragmentation 
function [24]: 

D"P(a;; 7) = (1 + 7)(2 + 7)(1 - x)x\ (13) 

and fit the free parameter 7 to the data. As in [16-18], we consider the data in the range 
0.18 < Xb < 0.94, to avoid the regions xb — and — ^ 1, where the calculation is 
unreliable. In fact, the parton- and hadron-level spectra presented in [13,16-18] become 
negative at very small and very large x, which is due to the presence of terms that have 
not been resummed yet, and to non-perturbative corrections, relevant especially at large 
X. When doing the fit, we set the scales in Eqs. © and (|7j) and the parameters entering 
the perturbative calculation to the following values: fi = fip = mz, yUo = A^of = ^b, 
mz = 91.188 GeV, rrih = 5 GeV, A^ = 200 GeV. In the considered range, we obtain: 
7 = 17.178 ± 0.303, with xVdof = 46.2/53 from the fit. 

If we try to compare the data with the considered Monte Carlo generators, we find 
that the default parametrizations of the hadronization models in HERWIG and PYTHIA 
are not able to reproduce the xb data. The per degree of freedom are, in fact, 
xVdof = 739.4/61 for HERWIG and xVdof = 467.9/61 for PYTHIA. A fit of HERWIG 
to a number of observables in e~^e~ annihilation was performed in [43]. Employing the 
parameters suggested in [43] gives x^/dof = 391.9/61 for the comparison with the xb 
data. 

In this paper we wish to reconsider the tuning of HERWIG and PYTHIA to the 
6-fragmentation data. We point out that we are aware of the fact that care must be 
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taken when changing the parameters of a Monte Carlo event generator, since this may 
spoil the agreement with the data sets that were used in the default tuning. The fits 
that we shall perform should not be seen as an attempt to redefine the parameters of 
HERWIG or PYTHIA; we would just like to investigate whether it is possible to change 
few parameters to describe the xb data better. Of course, it will be interesting to 
use the parametrization which we shall suggest to check how HERWIG and PYTHIA 
fare with respect to other observables. When doing the tuning, we concentrate on the 
parameters that determine the hadronization and leave unchanged the ones related to 
the perturbative phase of the showers and the quark masses. 

In HERWIG, we modify the values of CLSMR(l) and CLSMR(2), controlling the 
Gaussian smearing of the hadron direction with respect to the original constituent 
quarks; PLSPLT(2), which determines the mass distribution of 6-flavoured cluster de- 
cays; DECWT, affecting the relative weight of decuplet and octet baryons; and CLPOW, 
to which the heavy-cluster yield and the baryon/meson ratio are sensitive. The fitted 
values are: CLSMR(l) = 0.4 (default 0), CLSMR(2) = 0.3 (0), PSPLT(2) = 0.33 (1), 
DECWT = 0.7 (1), CLPOW = 2.1 (2). The values of CLSMR(l), PSPLT(2) and 
DECWT are indeed the same as in Ref. [43]. After the tuning, the agreement with the 
data is still not very good, but it is much better than with the default parametrization. 
We find xVdof = 222.4/61. 

As for PYTHIA, we change the values of the fragmentation parameters PARJ(41) 
and PARJ(42), which control the a and b parameters of the Lund symmetric fragmen- 
tation function 

hiz) oc -(1 - z)"exp(-6m2./^). (14) 

We also tune PARJ(46), which modifies the endpoint of the Lund function according to 
the Bowler hadronization model [46] ^: 

/^(^) °^ ;TTK (1 - ^)"exp(-&m|/2;), (15) 

where rriq is the quark mass. Our tuning gives the following values: PARJ(41) = 0.85 
(default value 0.3), PARJ(42) = 1.03 (0.58), PARJ(46) = 0.85 (1). After our tuning, 
PYTHIA matches the e^e~ data well, and we obtain x^/dof = 45.7/61 from the fit. In 
Tabled we summarize the parameters of HERWIG and PYTHIA that we have changed 
in our tuning. We have checked that our tuning works well also for the new model 
implemented in PYTHIA 6.3, as long as one runs it using options and parameters as 
they are defined in the new scenario (model 1) [47]. Running PYTHIA 6.3, we find 
X^/dof = 46.0/61 for the comparison with the xb data. In Figs. EJIHl we compare 
LEP and SLD data with HERWIG and PYTHIA respectively, using the default and 
tuned parameters. In both figures, we also show the results given by the resummed 
calculations of Ref. [13], convoluted with the Kartvelishvili model (solid), and denoted 

^For a quark of energy E moving on the z axis, the quantity z appearing in Eq. H14|l is the fraction 
oi E -\-pz taken by the hadron, while rriT is transverse mass, i.e. = m?' + '^Py- 

^See also the PYTHIA manual [20] for details on the implementation of Eqs. ((TH) and ifT^ . 
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Table 1: Parameters of HERWIG and PYTHIA hadronization models that we have 
tuned to improve the agreement with e^e~ data, along with the per degree of freedom. 



HKRWTfi 

J.XJ_;XVjVV -LVJr 


PYTHTA 


CLSMR(l) = 0.4 




CLSMR(2) = 0.3 


PARJ(41) = 0.85 


DECWT = 0.7 


PARJ(42) = 1.03 


CLPOW = 2.1 


PARJ(46) = 0.85 


PSPLT(2) = 0.33 




xVdof = 222.4/61 


xVdof = 45.7/61 



by 'NLO+NLL+Kart.' For the sake of comparison, we present the PYTHIA spectrum 
obtained using the parameters of Table ^ but without the rejection of non-angular- 
ordered showers. 

The xb distributions given by HERWIG and PYTHIA with the default parameters 
deviate significantly from the data points, as suggested by the values. The NLO-I-NLL 
calculation and PYTHIA after the tuning describe the data quite well. If we turn 
angular ordering off, PYTHIA is obviously not capable of reproducing the data. In 
fact, when tuning PYTHIA, we have used a version that vetoes non-angular-ordered- 
showers; if we allow non- angular-ordered branchings in PYTHIA, we should in principle 
even reconsider the tuning. Besides the comparison with the data, it is nonetheless 
interesting to observe the effect of possible non-angular-ordered showers: more events in 
the middle-low range 0.2 <xb< 0.6 and less events around the peak. In the following, 
we shall always run PYTHIA with forced angular ordering. 

The distribution yielded by HERWIG is instead somewhat broader than the data 
set. Even after the tuning, HERWIG is above the data at small and large xs and 
below the data in the neighbourhood of the peak. In any case, the comparison is 
greatly improved with respect to the default parametrization, and in principle even the 
HERWIG parameters which control perturbative h production can be further tuned. 
Furthermore, major improvements in the treatment of h fragmentation are present in 
HERWIG-|--|- [44], the new object-oriented C-I--I- version, which uses new showering 
variables and includes quark masses in the splitting functions [45]. As discussed in [44], 
one can obtain a quite good description of the SLD xb data fitting only the shower cutoff, 
corresponding to Qq in Eq. (jSI), and without any additional tuning of the hadronization 
model. The use of HERWIG-I--I-, whose first version simulates e^e~ annihilation, is 
however beyond the scope of this paper. 

Relying on the universality of the hadronization mechanism, we can use the parametriza- 
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Figure 5: Data from LEP and SLD experiments, compared with the NLO+NLL calcu- 
lation convoluted with the Kartvelishvili model (solid) and HERWIG 6.506, using the 
default parametrization (dashed) and our tuning (dotted). 




Figure 6: As in Fig. but comparing data and the NLO+NLL calculation with default 
(dashed) and tuned (dotted) PYTHIA 6.220. Also shown is the PYTHIA prediction, 
using our tuning of the hadronization model, but without rejecting the showers which 
do not fulfil angular ordering (dot-dashed). 

tions quoted in Table to predict the Xb distribution in H —>■ bb and t bW. In Higgs 
decay we shall investigate the cases rriH = 120 and 500 GeV; in top decay we shall as- 
sume rrit = 175 GeV and mw = 80.425 GeV. The other quantities will be set consistently 
with the values employed in the fit to the e^e~ data. 

In Fig. we plot the 5-hadron spectrum in H ^ bb processes, using HERWIG, 



13 



PYTHIA, and the NLL-resummed calculation in [18], convoluted with the Kartvelishvili 
model fitted as above. We have set the Higgs mass to mn = 120 GeV. PYTHIA fares 
rather well with respect the NLO+NLL calculation: although small discrepancies are 
present at very small and large xb and around the peak, the overall agreement looks 
acceptable. The behaviour of the curves given by HERWIG reflects what we found in the 
comparison with the e~^e~ data: even in Higgs decay, the xb spectrum is broader, hes 
above the NLO+NLL computation at large and intermediate xb, and below around the 
peak. Since matrix-element corrections to e^e^ bb were included when we did the fit, 
in principle the tuned parameters must be used only when matrix-element corrections 
to H ^ bb are turned on. However, for the sake of comparison, we show in Fig. [7| 
the Xb spectrum given by HERWIG without such corrections as well. The effect of 
matrix-element corrections is indeed what one should expect: we have more events at 
small Xb, say xb < 0.6, and less at large xb- In fact, as shown in Figs. ElandlH the 
implementation of hard and large-angle gluon radiation enhances the probability to have 
b quarks with smaller Xb, since the gluon is allowed to take a larger energy fraction. It 
is therefore reasonable that this is reflected by having more hadrons at small as well. 

We would like to consider also the case of a Higgs with a very large mass, in particu- 
lar rriH = 500 GeV. This example does not have a relevant application in the Standard 
Model, where the branching ratio of H —>■ bb is tiny for such a high-mass value, but 
it could become more important in models with an extended Higgs sector. Another 
interesting aspect is to investigate whether our predictions are still in reasonable agree- 
ment for masses well above the Z-boson mass. For ttih = 500 GeV, width effects may 
become important; however, for the sake of comparison with the resummed calculation, 
which does not implement the Higgs width, we run HERWIG and PYTHIA for Th = 0. 
The spectra shown in Fig. |H1 look broader than the 120 GeV case, but the comparison 
is rather similar: PYTHIA yields the highest spectrum at small xb and HERWIG at 
intermediate and very large xb- The effects of matrix-element corrections to H ^ bb 
in HERWIG has an impact that is similar to the one already observed if Fig. [7| more 
events at small xb and less around the peak. We have finally checked that the impact 
of the inclusion of finite-width effects in HERWIG and PYTHIA is of the order of 1%. 

In Fig. IHl we compare the i?-hadron spectrum in top-quark decay, using HERWIG, 
PYTHIA and the perturbative NLO+NLL resummed calculation performed in [16, 17], 
convoluted with the Kartvelishvili non-perturbative model. PYTHIA is able to repro- 
duce rather well the peak given by the resummed calculation, while its spectrum lies 
below it at xb ^ 0.7 and above it at > 0.95. The spectrum given by HERWIG looks 
instead slightly different: it is below the resummed calculation in most of the xb range, 
i.e. Xb ^ 0.9, and above at very large xb- 

In [18] the B spectra in Z, top and Higgs decays were compared using NLL resummed 
calculations. In Figs. ^1 and ^2 we show the same comparison, but running HERWIG 
and PYTHIA; in order to investigate the effect of the particle spins, we also show the 
result in e^e~ annihilation, but for ^/s = 120 GeV, the chosen value for the Higgs mass. 
Figs. ITUlandfTTlexhibit features similar to those of Ref. [18]: the B energy distribution in 
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Figure 7: i?-hadron spectra in H ^ bb processes, according to a NLO+NLL calcula- 
tion with the Kartvelishvili non-perturbative model (solid line), PYTHIA (dotted) and 
HERWIG with (dot-dashed) and without (dashed) matrix-element corrections. We have 
set niH = 120 GeV and used the tuning in Table [TJ 
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Figure 8: As in Fig. [3 but for a Higgs mass mu = 500 GeV, and with no width effects 
in HERWIG and PYTHIA. 



top decay is shifted toward large xb values; in Higgs decay we have more events at small 
xb] the spectrum from Z decay lies between the two. Setting mz = rriH, as observed 
in [18], makes the spectra from H and Z decays more similar, but small discrepancies 
are still present. The comparison between top and Higgs decays is indeed quite relevant, 
since the different B spectra may help to distinguish between 6-flavoured hadrons in the 
process pp tiH, followed hj H ^ bb, one of the most important channels for Higgs 
studies at the LHC [7]. 
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Figure 9: i?-hadron spectra in top decay, for rrit = 175 GeV, according to a NLO+NLL 
computation convoluted with the Kartvehshvih model, (solid line), HERWIG (dashed) 
and PYTHIA (dotted). 




Figure 10: Comparison of i?-hadron spectra m Z hh for uiz = 91.188 GeV (solid) 
and mz = 120 GeV (dot-dashed), for top decay (dashed) and for H ^ bb (dotted), with 
rriH = 120 GeV. Results given by HERWIG, with matrix-element corrections included 
in all plots. 
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Figure 11: As in Fig. [El but using PYTHIA. 



6 Results in moment space 

In this section we wish to present resuhs in moment space, where the moments F^v of a 
function 1/F dT/dz are defined as follows: 



In the case of the histograms given by a Monte Carlo program, the integral in ()16p will 
be approximated to a discrete sum over the bins. In Ref. [48], the DELPHI collaboration 
presented the moments for B production in e'^e~ annihilation. From the point of view 
of resummed calculations, working in moment space [26,49] presents several advantages. 
In A^-space, convolutions become ordinary products, and the relation between parton- 
and hadron-level cross sections becomes: 



where cr^ and are the moments of the b and B cross sections, and D'^ the A^-space 
counterpart of the non-perturbative fragmentation function. Therefore, there is no need 
to assume any functional form for the non-perturbative fragmentation function in x 
space. Moreover, resummed calculations are well defined in A^-space, and we do not 
have the problem exhibited by the x^-space spectra, which become negative at small or 
large xb- 

In [17,18], the parton- and hadron-level moments for t —* bW and H hb processes 
were presented, using a resummed NLL perturbative calculation. In this paper, we wish 
to compare experimental data and predictions, employing HERWIG and PYTHIA as 
well, tuned as in Table [T] 




(16) 





(17) 
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In Table El we quote the data from DELPHI, the moments yielded by HERWIG and 
PYTHIA in Z, t and H decays, and, for the sake of comparison, the results, already 
published in [17, 18], yielded by NLL-resummed calculations. As for the NLL computa- 
tion, we present two sets of results: we multiply the moments of the perturbative cross 
sections (widths) by the moments of the Kartvelishvili model, denoted by K]\f in Ta- 
ble 121 fitted to LEP and SLD x^-space data, and by -D^, extracted from the DELPHI 
iV-space data [48]. The moments yielded by HERWIG and PYTHIA in e~^e~ anni- 
hilation are consistent, within the error ranges, with the ones measured by DELPHI. 
Although problems are present when fitting the xb data from LEP and SLD, it is re- 
markable that HERWIG is compatible with the DELPHI moments within one standard 
deviation. Moreover, we observe that the moments given by the NLL calculation are 
rather different according to whether the fit is made directly to the A^-space DELPHI 
data, or by fitting in xb space and then calculating the moments. In fact, the x^-space 
fit was performed in the range 0.18 < xb < 0.94, while the moments are computed 
integrating over the full range < < 1. Neglecting data at small and large xb in the 
fit may cause an incorrect determination of the moments, which do depend on the tails 
of the distributions [26,49]. 

The results for top and Higgs decay exhibit similar features to the xb spectra. In 
top decay, PYTHIA is very close to the NLL calculation which uses D'^f extracted from 
the DELPHI data, while HERWIG, whose predictions are shifted toward larger xb, 
gives larger moments. For H ^ bb and mH = 120 GeV, PYTHIA and HERWIG, after 
matrix-element corrections, give moments which are compatible within 1%, and are 
closer to the NLL result which uses the Kartvelishvili model rather than the moments 
Z)]^. For niH = 500 GeV, HERWIG looks somewhat closer to the NLL result, while 
PYTHIA gives moments which are a few per cent smaller. Matrix-element corrections 
to H ^ bb in HERWIG have the effect to reduce the values of the moments of about 
2-5% for both rriH = 120 and 500 GeV, which is understandable, since more small-xs 
events are generated after the inclusions of such corrections. 



7 Conclusions 

We have studied bottom-quark fragmentation in e'^e' annihilation, as well as in top 
and Higgs decay, using NLL-resummed calculations based on the fragmentation function 
formalism and Monte Carlo event generators, such as HERWIG and PYTHIA. Monte 
Carlo parton showers simulate multiple radiation in the soft or coUinear approximation, 
and need to be provided with matrix-element corrections to simulate hard and large- 
angle parton radiation. While PYTHIA contained the corrections to all three considered 
processes, we had to implement matrix-element corrections to H —>■ bb in HERWIG, since 
they are not yet present in HERWIG. At parton level, the effect of such corrections is 
that there are more b quarks simulated with a smaller energy fraction. 

We have fitted HERWIG, PYTHIA and the Kartvelishvili hadronization model, the 
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Table 2: Moments from DELPHI [48], and moments in e+e annihilation, Higgs (H) 
and top (i) decay, using NLL resummed calculations, HERWIG (HW) and PYTHIA (PY). 
In if — > 66 , we consider the Higgs mass values tjih = 120 GeV and 500 GeV, and quote 
HERWIG results with (c.) and without (no c.) matrix-element corrections. The NLL results 
are obtained by extracting D^^ from the DELPHI data, and using the moments of the 
Kartvelishvili model, fitted to x_b data, as in Figures El and El 





(x) 




(a;3) 




e+e data 


0.7153±0.0052 


0.5401±0.0064 


0.4236±0.0065 


0.3406±0.0064 


e+e- NLL 


0.7801 


0.6436 


0.5479 


0.4755 




0.9169 


0.8392 


0.7731 


0.7163 


e+e~ = a^K^ 


0.7027 


0.5251 


0.4066 


0.3225 


e+e- HW 


0.7113 


0.5354 


0.4181 


0.3353 


e+e" PY 


0.7162 


0.5412 


0.4237 


0.3400 


if-dec.(120) NLL r% 


0.7580 


0.6166 


0.5197 


0.4477 


//-dec. (120) = r%D''^ 


0.6950 


0.5175 


0.4018 


0.3207 


i?-dec.(120) = T^Kn 


0.6829 


0.5030 


0.3858 


0.3036 


/7-dec.(120) HW c. 


0.6842 


0.5036 


0.3877 


0.3076 


//-dec. (120) HW no c. 


0.6961 


0.5177 


0.4011 


0.3197 


i/-dec.(120) PY F^ 


0.6876 


0.5080 


0.3913 


0.3099 


i/-dec.(500) NLL F^^ 


0.6858 


0.5265 


0.4255 


0.3545 


//-dec. (500) F^ = T%D''^ 


0.6288 


0.4418 


0.3290 


0.2539 


//-dec. (500) F^ = T^Kn 


0.6178 


0.4295 


0.3158 


0.2404 


//-dec. (500) HW c. F^ 


0.6184 


0.4279 


0.3146 


0.2406 


//-dec. (500) HW no c. F^ 


0.6286 


0.4389 


0.3245 


0.2491 


//-dec. (500) PY F^ 


0.6044 


0.4152 


0.3036 


0.2307 


t-dec. NLL F^ 


0.7883 


0.6615 


0.5735 


0.5071 


t-dec. NLL F^ = T^D"}^ 


0.7228 


0.5551 


0.4434 


0.3632 


t-dec. F^ = T%Kn 


0.7102 


0.5397 


0.4257 


0.3439 


t-dec. HW F^ 


0.7325 


0.5703 


0.4606 


0.3814 


t-dec. PY F^ 


0.7225 


0.5588 


0.4486 


0.3688 



latter used in conjunction with the NLL-resummed calculation, to e+e data on the B- 
hadron energy fraction xb, from ALEPH, OPAL and SLD. We found that the Kartvel- 
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ishvili and the PYTHIA string model are able to fit the data quite well, while the 
HERWIG cluster model is only marginally consistent. The tuning of the Monte Carlo 
programs turned out to be essential, since the default parametrizations fare rather badly 
against the measured xb spectrum. 

Relying on the process- independence of the hadronization, we have predicted the xb 
spectrum in Higgs and top decays, using the tuned models, and found results which 
reflect the quality of the flt to the e+e^ data. PYTHIA and rcsummcd calculations 
give consistent results, while the HERWIG predictions are broader and compatible only 
in some xb ranges. The effect of matrix-element corrections to if 66 is reasonable, 
since more events are generated via the exact amplitude at small xb, which corresponds, 
at parton level, to hard and large-angle gluon radiation. Comparing the spectra of B 
hadrons in Z, H and top decays has given similar results to the ones already found in 
the framework of perturbative fragmentation functions: the -B's from top quarks are the 
hardest, the ones from H are the softest, and the ones from Z lie within the two. 

Finally, we have considered data on the moments of the B cross section measured 
at LEP by the DELPHI collaboration. We have extracted the moments of the non- 
perturbative fragmentation function, and compared the A^-space cross sections given by 
the NLL computation to those yielded by HERWIG and PYTHIA. We found that tuned 
HERWIG and PYTHIA are able to reproduce the DELPHI moments, within the quoted 
error range, even though HERWIG had problems with fltting the e^e~ Xb data. The A^- 
spacc predictions in top decay reflect the behaviour of the x^-space results: PYTHIA is 
closer to the NLL calculation, and HERWIG yields larger moments. In if — > 66, all given 
moments are pretty similar ior uih — 120 GeV, while, for mn — 500 GeV, HERWIG 
is closer to the NLL prediction and PYTHIA yields smaller results. Matrix- element 
corrections to HERWIG simulations of ii ^ 66 decrease the moment values, which is 
due to the smaller fraction of events at large xb- We have also pointed out that the 
moments based on an NLL x^-space flt, making use of the Kartvelishvili fragmentation 
model in a range which discards very large and very small values oi xb, are different 
from the ones obtained after fitting the moments experimentally measured by DELPHI 
directly. 

In conclusion, we believe that our studies and tunings can be a useful starting point 
to address 6-quark fragmentation at the Tevatron and LHC, especially for the analyses of 
top-quark properties, such as its mass, and for the Higgs studies in the H ^ bb channel. 
It will also be interesting to use the parametrizations which we have proposed to study 
other hadron-level observables not necessarily related to 6-quark fragmentation. More- 
over, while we deliberately did not aim at fitting any perturbative parameter entering in 
the Monte Garlo simulation or in the resummed calculation, it will be nonetheless cum- 
bersome comparing parton-level observables and investigating whether further tuning 
may improve the comparison with the e"'"e~ data and affect the hadron-level predictions. 
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